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SELF -STARTING MULTISTEP METHODS FOR THE NUMERICAL 


INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS 

By William A, Mersman 
Ames Research Center 


SUMMARY 


Classical, multistep; predict or- correct or procedures for the numerical 
integration of systems of ordinary differential equations are generalized to 
provide compatible; self -starting methods. Explicit algorithms and tables of 
numerical coefficients are presented. 


INTRODUCTION 


The numerical integration of systems of ordinary differential equations 
on modern automatic computers is usually accomplished by means of so-called 
multistep methods; particularly the predictor-corrector methods associated 
with the names Adams; Bashforth; Moulton; Stormer; and Cowell, It is usually 
assumed that these methods are not self- starting; and recourse is had to 
single-step methods like that of Runge-Kutta to obtain starting values. This 
leads to cumbersome computer programs requiring what amounts to unessential 
tallying to determine whether enough starting values have been obtained. 

The purpose of the present report is to derive simple generalizations of 
the classical predictor-corrector formulas that immediately yield compatible 
self -starting procedures that produce all the required backward differences 
directly from the initial conditions. 


STATEMENT OF THE PROBLEM 


The problem is to devise a self -starting; multistep procedure for the 
numerical solution of the initial value problem 


dx 

dt 





x(t Q ) = x 0 

y(t Q ) = y Q 


(i) 


(2) 



at the discrete, equally spaced points t n , n = 1, 2, 3, • . . The vari- 
ables x, y, and f are vectors, all of the same (finite) dimension. 


that 


Let the common interval of the independent variable be denoted by 

t n - ^0 + nh 


h, so 


and introduce the usual notation 

x n = x(t n ) 

yn = y(t n ) 

f n = f(x n ,y n ,t n ) 

The index, n, will be restricted to integral values, usually positive, 
although negative values will be introduced in some of the starting procedures 
to be discussed later. 


The general problem, then, is to devise algorithms for calculating x n , 
y n , f n , for n = 1, 2, 3, . . . , given simply the differential equations (l) 
and the initial values ( 2 ) . The theory for first-order systems is obtained by 
ignoring the variable, x, throughout the general theory. 

The procedure to be used is the conventional one of approximating the 
function, f, by a polynomial in t of degree q. The problem is then split 
into two: the forward integration problem and the starting problem. 


The Forward Integration Problem 

The problem of integrating forward one step will be solved by means of 
backward difference formulas of the Adams type. Here it is assumed that t n , 
x n , yn> fn> an< 3- the first q backward differences of f n : 

'N 

V°fn = fn 
Vfn = fn “ f n-i 

, . . > (3 

= V^ _1 f n " V^ -1 fn- 1 

k = 1(1) q. 

are known. Several algorithms will be derived for computing all these 
quantities at t = t n+1 . 


The Starting Problem 

Before the forward integration algorithms can be applied, it is necessary 
to compute initial values of the backward differences of f at some point, 
preferably at t = t 0 . This will be done by developing iterative algorithms 
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for the backward calculation of x_ n , y_ n , f-n n = l(l)q, from which the 

backward differences of f Q are easily calculated. 

It sometimes happens that the initial values lie near a singularity. In 
this case the backward starter may fail. For this reason, iterative forward 
starters will also be derived for the calculation of x n , y n , f n , at 
n = l(l)q. 

In either case the final "output" of the starting procedure will be x Q , 
y G , f o j and the first q backward differences of f Q , so that the starter 
will be compatible with the forward integration procedures. 

In the following section the basic backward difference equations will be 
derived. These are generalizations of the equations usually ascribed to 
Adams, Bashforth, Moulton, Stormer, and Cowell (ref. 1, chs. 5 and 6). It is 
the generalization of the classical formulas that makes self -starting proce- 
dures possible. 


GENERALIZED BACKWARD DIFFERENCE EQUATIONS 


Normal Form 

The basic difference equations, relating Vy, Vx, and V 2 x at t = t j to 
f and its backward differences at t = t n , both j and n being arbitrary, 
are 

oo 


Wj - h 7n- j 

k=o 
00 

^ ff n-j,k vkf n 
k=o 

CO 

VXj = hy n + h 2 ^ (a n _j^k+i - 7o,k+l)v k fn 


V 2 Xj = h 2 


(10 

(5) 

( 6 ) 


k=o 


These are written formally as infinite series to simplify certain index manip- 
ulations. They terminate and are exact whenever f is a polynomial in t. 

The proof is a straightforward generalization of Henrici’s (ref. 1, 
pp. 191-19^ and 290-293)* Write yj as a Taylor's series with remainder 
centered at t j_i (ref. 2, p. 95 ) • 


y j = y j-i +h + rh)dr 


Now note that 


ti-i + Th = tn - (n - j + 1 - T)h 
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Approximating f 
pp. 190-191) 


by means of Newton's interpolating polynomial (ref. 1, 


00 



k=o 


•where the symbol 



r(s ± 1) 

kir(s + 1 - k) 


yields equation (4), with 7 given by 


7 P,k 


(-D* 


rr" 



(7) 


Equation (5) is obtained similarly by writing the Taylor's series for 
Xj and Xj_ 2 centered at tj_ x 


x .i = x .i-i + 


hyj-x + h 2 r (l - t ) f (t + T h)dT 

Jo 


X j-2 


X j 


-1 



T)f(tj_l - rh)dr 


adding and inserting Newton's interpolating polynomial yields equation (5) , 
with a given by 


°p,k 




+ 1 
k 




+ 1 + t\ 

k / 


dT 


( 8 ) 


Before deriving equation (6) it is convenient to discuss equations (4) 
and (5) and some of the properties of 7 and cr. 

Equation (4) with j = n + 1 is the Adams -Bashforth predictor (ref. 1, 
pp. 192-193)* Henrici uses the notation 

7 k = 7-x, k 

Equation (4) with j = n is the Adams-Moulton corrector (ref. 1, 
pp. 194-195) • Henrici uses the notation 

7k* = 7o,k 

Equation (5) with j = n + 1 is the Stormer predictor (ref. 1, pp. 291- 
292 ). Henrici uses the notation 


Equation ( 5 ) with j = n is the Cowell corrector (ref. 1 , pp. 292-293). 
Henrici uses the notation 

a k* = a o,k 

As will be seen later , equations (4) to ( 6 ) with j = n - q(i) n - 1 
form the basis for the self -starting procedures to be developed. 


The most important property of 7 and a is that each row is the first 
backward difference of the preceding row: 


7p+i,k - 7p,k - 7p,k-i 

°p+i,k = CT p,k " CT P^k-i 


k > 1 


(9) 


while 7 p 0 = cTp^o - 1# These follow immediately from the definitions , equa- 
tions ( 7 ) * (8), and the well-known recurrence relation for the binomial coef- 
ficients : 



Equations ( 9 ) can be rearranged in the useful form 


v j7 n -j,k+i “ 7 n -j,k 

v j CT n-j,k+i = G n-j,k 


( 10 ) 


where the subscript on Vj is used to emphasize that V here operates on j, 
not on n or k. 


The tables of 7 and o at the end of the report were computed by means 
of the definitions, equations ( 7 ) and ( 8 ), for the first row (p= -l), and the 
difference equations ( 9 ) for subsequent rows. 


To return to the derivation of equation ( 6 ), note first that x bears 
the same relation to y as y does to f. Hence, we can write equation (4) 
in the transliterated form 

00 

m=o 

Writing equation (4) with j = n gives 


00 



r=o 

and taking the (m - l)st backward difference gives 

00 

- 'Yto,****- 1 ** 

r^o 


v*y 


n 
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Substituting in equation ( 4 a) and rearranging gives 

oo k 


Vn + h2 ^ ^ 7 0 ,r7n-j,k+i-r v ^n 


( 11 ) 


k=o r=o 


Applying the backward difference operator V- and using equations ( 10 ) gives 

J 


oo K 

V2 Xj =h 2 ^ ^7 0 , r 7 n _j, k _ r V k f n 


k=o n=o 

Comparison with equation (5) gives the important identity 

k 


Hence, 


k 


c p,k ~ ! ^ 7 p,r 7 p,k-r 
r=o 


7 o,r 7 p,k+i-r - a p,k+i “ 7 o,k+i 


( 12 ) 


r=o 


and this reduces equation (ll) to equation (6)* Q.E.D. 

Equation (l2) constitutes a valuable check on the tables of 7 and o and 
has been used. In addition, for small values of k, it is convenient to have 
the explicit formulas 


7 p,o 

= 1 

7 P,i 

= -P " 

7 P,2 

C\J 

ft 

H| CVJ 
II 

a p,o 

= 1 

CT P, 1 

= -p - 


_ p(p • 

G p,2 

2 


1 

2 


12 


(13) 


12 


m 


Equation (6), with j = n + 1 or n, constitutes a new predictor- 
corrector scheme, to be discussed later, which appears to have some advantages 
over the St ormer -Cowell scheme of equation ( 5 )* Eor j = n - q(l)n - 1, equa- 
tion (6), like equations ( 4 ) and (5), forms the basis for an iterative 
starter. 

This completes the derivation of the basic backward difference equations 
in the normal form. Closely related to these are similar equations using the 
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first and second sums. Astronomers have long used these concepts (ref, 3 ), 
and further confirmation of their effectiveness is offered by Henrici (ref. 1, 
PP* 327-339) • These summed forms will be derived in the next section. 


Summed Form 


The basic difference equations, in the summed form, relating x and y at 
t = tj to f, its backward differences and its backward sums at t = t n , both 
j and n being arbitrary, are 


00 

yj = An + ^ ^ ?n- 
k=o 

00 

x j = B n + (J - n - l)hA n + h 2 ^ 


(15) 


k=o 

where the backward sums A^. and B n are defined by these same equations with 
j = n: 


A n - yn - 

k=o 

00 

B n = x n + hA n - h 2 ^ Oo^k+a^n 


(16) 


k=o 

The names first and second sum for A n and B n , respectively, are justified by 
the property 


VA n = hf n 
VB n = hA n 


(17) 


The proof is quite straightf orward. Replace the coefficient, y, in equa- 
tion (4) by its equivalent value from equations (10 ) : 


✓ 00 \ 
V j y j = V J( h 7 y n-j,k+iV k fn 


k=o 

where, again, the subscript j is affixed to V to emphasize that the opera- 
tor, especially in the right member, acts on j rather than n or k. This is 
a simple, linear difference equation whose solution is obviously the first of 
equations ( 15 ), in which A n is an arbitrary constant vector, independent of 
j, whose value is determined by setting j = n, equations (l6) . 

If y^ is now eliminated from equation (6) by means of the first of 
equations (l 6), and if the coefficient, a, is replaced by means of equa- 
tions (10), the result is again a simple, linear difference equation: 
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r 


v j x j 


= hA n + ( h 2 ^ cr n _ 

V k=o 


Jl “ [_> w n-j,k+2 


V^n 


whose solution is clearly the second of equations ( 15) , where Bn is another 
arbitrary constant vector independent of j, whose value is determined by 
setting j = n (eqs. (l6)). 


Equations (17) follow immediately from equations (l 6) on taking the back- 
ward difference and applying equations ( 4 ) and (6) with j = n. 

As will be seen in later sections, equations (l 6) will be used only once, 
in connection with the starting procedure, to obtain initial values, A 0 and 
Bo, of the first and second sums. Subsequently, equations (15) with j = n+1 
will yield the summed versions of the Adams-Bashforth and Stormer predictors, 
while equations (15) and (17), with n replaced by n+1 and j = n + 1, will 
yield the summed versions of the Adams-Moulton and Cowell correctors. 


Before discussing starting procedures, it is necessary to convert the 
basic difference equations ( 4 ), (5), (6), ( 15 )* and (l 6) to the backward ordi- 
nate form, in which earlier values of f are used rather than its backward 
differences . 


GENERALIZED BACKWARD ORDINATE EQUATIONS 


The backward difference equations of previous sections have been written 
formally as infinite series merely for manipulative convenience. In practice, 
it is assumed that f is approximated by a polynomial in t of degree q, 
say, and the series are all truncated at k = q. If, then, the backward dif- 
ferences in equations (15) and (l6) are eliminated by means of the well-known 
identity (ref. 1, p. 190) ^ 

^fn =^(-D m (m) f n-m 
m=o 


and the order of summation is reversed in the resulting double sums, the 
following backward ordinate equations are obtained in summed form: 


1 

V 


yj An + h 7 ^ 7q,n- j,m^n-m 
m=o 

x j = ®n + " n ~ + k 2 X^n-^n- 


m=o 


where the coefficients 7 and a are 


(18) 
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7q,p,m ~ (“l) 7p,k+i 



m 

k=m 

= ( -1 ) °p,k+2 

k=m 

Again the sums A n and Bn are given by setting j = n: 

a 

A n = Yu ^ I ^q, o,m^n-m 

m=o 


( 19 ) 


a 

B n = x n + - h 2 ^tfq.,o,m f n-in 


m=o 


( 20 ) 


Eliminating A n and B n from equations (l8) and (20) then gives the normal 
form of the backward ordinate equations : 

CL 


j - y n - h^T 


a q,n- j^m?n-m 


m=o 


x 


3 = x n - (n - j)hy n + h 2 ^ 


Pq.,n- j^m^n-m 


m=o 


( 21 ) 


where the coefficients a and 0 are 


a q.,p,m - 7q_, o,m " 7q.,p,m 

Pq,p,m = a q.,p,m “ a q.,o,m + p7 q.,o,m . 


( 22 ) 


The tables of y 9 a, a, and [3 at the end of the report were computed 
using equations ( 19 ) and (22) together with the easily proved relations: 


7 <i,p,° 

tfq.,p,o 

7 q.,p,q. 

a q.,p,q. 

7 q+i,p,m 

a q. +1 >p>m 


7 P-I,q.+1 " 1 


Q'p-l^q +2 + p + 1 


(”l)^ 7 Pjq+i 


( - - 1 -) < ^ a p,q.+2 


7 + (-l) m + lN l 

7 q.,p,m + L -l; ^ m J 

? P,a-te 

+ (-l) m ( q m ^ 

1 a p*q+ 3 
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Taking f(t) = t r in equations (21 ) , r = 0(l)q, gives the identities 

^m r a,^ p#m = % — - , ^ 7 q.,p,m = 7 P,i 

m=o m=o 

q. q. 

" (r + l)(r + 2) ’ ‘ <>p,a 

m=o m=o 

which were used as a final check on the tables. 

All the basic formulas, summed and normal, in backward difference and 
backward ordinate forms, have now been derived. The remaining sections of the 
report are devoted to the presentation of general and specific algorithms for 
starting and continuing the integration. 


BACKWARD STARTING MjGORTEHMS 


A backward starting algorithm is one that produces the backward differ- 
ences of f 0 , up to order q, given merely the initial conditions, equa- 
tions (2), and, of course, the differential equations (l). This is equivalent 
to an algorithm that will produce backward values of the ordinates, f^p, for 
p = l(l) q; it is then a simple matter to obtain the backward differences at 
f Q . (See the appendix for details.) 


The desired algorithm, in normal form, is implicitly contained in equa- 
tions (21) with n = 0, j = -p: 


y-p = Yo ~ h 


m=o 


x_^ = x n - phy^ + h‘ ) Pq^p, m f-m 

m=o 

f_p = f (to - ph,x_ p ,y_ p ) 


(23) 


for p = l(l)q. This is a set of 3q implicit equations for the 3q unknowns 
y -P > x -p> f -p* 


If an approximate solution is known., an improved solution is readily- 
obtained by iterating equations (23); inserting "old" values in the right mem- 
bers produces "new" values in the left members. Collatz (ref. k, pp. 99 -IOI) 
exhibits the a for q = 2, 3 and. discusses the convergence of the iteration, 
but does not derive the general formulas. Two methods of obtaining an initial 
approximation will now be given. 
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Bootstrap Starter 


Since only one value of f, namely f Q , is known initially, the simplest 
possible procedure is to take q = 0 in equations ( 23 ). Setting p = 1 gives 
the "predicted" values of x, y, and f at t-i* Then setting q = 1, p = 1, 
gives corrected values. Keeping q = 1 and now setting p = 2, gives pre- 
dicted values at t- 2 . This bootstrapping procedure can be repeated until the 
desired value of q is reached. The explicit algorithm is 


Predictor 


p = 1 ( 1 ) 0 . 


-\ 


p-i 


y~p - k ^ a p-i,p,m-^-m 

m=o 

p-l 

x -p = x o - phy D + h2 

m=o 

f_p = f (t-p,x_p,y-p) 

Multiple corrector k = l(l)p 

P 

y-k “ ^o " ^ ^ a p,k,m^-m 
m=o 

P 

x -k ” x o " ^ Pp,k,m^-m 


m=o 


f_k = f(t-k,x-k,y-k) 


\ (24) 


J 


Experienced computer programmers will recognize this as a simple, nested DO 
LOOP. 


To make the procedure more tangible, the first few algorithms are written 
explicitly below: 


p = 1 ; predictor 


k = 1 ; corrector 


y-i = y 0 - hf 0 

h 2 

x-i = x 0 - hy Q + -- f o 

f-x = f (t-i,x- 1 ,y- 1 ) 

y-i = y 0 - | ( f o + f-i) 

x-1 = Xq - hy 0 + V ( 2f ° + f -i) 
o 

f-l = 
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p = 2; predictor 


y-2 = y 0 - 2hf-i 


x_2 = x 0 - 2hy 0 + ^ (2f 0 + 4f_x) 
f-2 = f(t- 2 ,x- a ,y_ a ) 


Multiple corrector 
k = 1 


y-x = y 0 -^(5fo + 8f-i - f-2) 

h 2 / 

X— 1 = Xo - hy 0 + ^ (7fo + 6f-x 
f-x = f(t-i,x-x,y_i) 


f- 2 ) 


k = 


2 


'y-2 = y 0 - | (f 0 + *tf-x + f-2) 

h 2 

X— 2 = Xq - 2hy 0 + ~ (2fo + 4f-x) 


[ f-2 = f(t-2,x-2,y- 2 ) 


p = 3; predictor y- 3 “ y 0 - 7* (3*o + 9 ?-z) 

X— 3 = Xo - 3hyo + ^ (9fo + lfif-1 + 9f- 2 ) 

f -3 = f(t-3,x_3,y_ 3 ) 


Multiple corrector 
k = 1 


' y-x = yo - ^ (9f 0 + 19f-x - 5f-2 + 

x -x - x 0 - hy Q + (97^0 + ll^f-i 


f-3) 


- 39f-2 


+ 8f- 3 ) 


l f_x = f(t-x,x-x,y-i) 


k = 2 


y-2 

X-2 

L f -2 


- yo - ^ (fo + ^-f-x + f-2) 

= x c - 2hy 0 + ^ (28fo + 66f-x 

45 

= f(t-2,x-2,y- a ) 


6f — 2 + 2f- 3 ) 


k = 3 


y -3 - yo - ^ (3fo + 9f-x + 9f-2 + 3f-3) 

x-3 = x Q - 3hy 0 + ^r~ (351fo + 972f-x + 243f- 2 + 54f- 3 ) 

360 " 

f-3 = f(t-3,x- 3 ,y_ 3 ) 


The bootstrap starter seems to be quite efficient in practice., but it is 
awkward and space- consuming when programmed for automatic computers, because 
of the multiplicity of matrices and algorithms required. This suggests the 
following logically sampler method. 
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Iterated Starter 


The bootstrap starter is essentially an efficient method of obtaining 
first approximations for use in the right members of equations (23), which are 
then iterated. A logically simpler, but less efficient method is to initial- 
ize by setting f_ m = f Q , m = l(l)q, in the right members, and then iterate 
the single set of equations (23) • 


Backward Starter, Summed Form 

Starting with equations (l 8 ) and ( 20 ), instead of (23) and again setting 
n = 0 , j = -p gives the implicit, summed form of the backward starter: 

p = 1(1) q. 


y-p A 0 + h ^ yq,p,ra^-m 

m=o 

QL 

x_p = B 0 " (p + l)hA 0 + h ^ ^q^p^m^*- 


m 


m=o 


f -p - f (" b -p^ x -p^y-p) 

q. 

A o = yo h ^^q,o,m f - 




(25) 


■m 


m=o 


B 0 = x 0 + hA 0 - h" 


'll. 


f 

nr -m 


j 


m=o 


Initial values can be obtained by the obvious bootstrap procedure or by start- 
ing with f_ m = f Q , m = l(l) q. 


A o = y Q + | f o 


ti 


B 0 = x Q + hA 0 - — f 0 = x Q + hy Q + ^ 
Equations (25) can then be iterated until they converge. 


-5A h 2 


h f , 


This algorithm is mentioned mainly for the sake of completeness. The 
principal reason for introducing the first and second sums is to obtain better 
control of the accumulated round-off error during a long integration, but this 
consideration may be irrelevent to a starting procedure. 

Before turning to the subject of forward starting procedures, it may be 
noticed that the backward starter (eqs. (23)) produces the ordinates f-p, 
p = l(l)q. These are easily converted to backward differences at f Q (see the 
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appendix) . If summed forms of the forward integration procedure are to be 
used subsequently, the initial values of A 0 and B 0 are easily computed from 
equations (l6) , with n = 0. Thus even here there is no compulsion to use the 
summed form of the starter. 


FORWARD STARTING ALGORITHMS 


The derivation of forward starting algorithms is almost trivial. Each 
of the backward starters discussed previously becomes a forward starter by 
means of the simple transf ormation 


h -h 
x -m x m 

y-m y m 

f -> f 
- 1 -m 

This produces the forward ordinates, 111 = 0(l)q. 

The conversion to backward differences at tq is straightforward. How- 
ever, to do this on an automatic computer would involve either losing all 
information at the points between t 0 and tq, or else increasing the storage 
requirements excessively. Furthermore, either choice leads to a procedure 
that is different, in its external appearance and mode of usage, from the 
backward starter. 

A preferable procedure is to compute the backward difference table at 
tq and then extend it back to t 0 by holding V^-f constant; this, of course, 
is consistent with the starting procedure. This is easily done, and program- 
ming details are given in the appendix. 


Choosing the latter alternative provides computer programmers with a 
battery of starting procedures, forward or backward, bootstrap or iterative, 
in normal or summed form, with identical external appearance. In every case 
the input data consist of the initial conditions, and the -output data consist 
of the table of backward differences (and sums) at the initial point, in a 
form compatible with the forward integration procedures to be discussed next. 


PREDICTOR-CORRECTOR ALGORITHMS FOR FORWARD INTEGRATION 


The purpose of this section is to present a variety of algorithms for the 
forward integration from t n to t n +j_* Specifically, it is assumed that the 
input consists of t n , x n , y n , f n , and the first q backward differences of 
f n , together with the sums A n and B n when appropriate. The output is to be 



the same list of quantities at t n +i. The combination of any of these algo- 
rithms with any of the starters provides the complete solution of the initial- 
value problem. 

All the algorithms to be presented involve the use of a predictor fol- 
lowed by a corrector, requiring two calculations of f (x,y,t) . Conflicting 
philosophies regarding the need for a corrector can be found in references 3 
and 5. The general consensus among automatic computer users seems to favor 
the use of one corrector. In the algorithms given below the user can, of 
course, omit the corrector if he chooses. 

The backward difference equations (4) through (6) and (15) through (17) 
give predictor formulas when j = n + 1. Taking j = n and then replacing 
n by n + 1 yields corrector formulas. In every case, if both predictor and 
corrector are truncated at the same value k = q, then subtracting the pre- 
dictor from the corrector and using the recurrence relations (eqs. (9 )) for 
the coefficients 7 and a gives a shorter formula for the corrector. 

Throughout, predicted values are indicated by an asterisk (*) • 


Normal Form 

First-order system .- The Adams -Bashforth predictor is 

q. 

jSm. = yn + h X 7 - l ’ k7kfn (26) 

k=o 

and the Adams-Moulton corrector is 

y„+i = 4+i + h ’'-i,a v ' 1+1 4 + i (27) 

Simple second- order system .- If the first derivative, y, does not occur 
in the differential equation and is not required, the formulas for x are: 

Stormer predictor ^ 

x n+i = 2x n - x n-i + h2 a -i,k vkf n 

k=o 

Cowell corrector * q+i * 

x n+i = x n+i + k °*-i^q^ -^n+i (29) 

General second-order system .- When the first derivation, y, is present, 
equations (4) and (6) yield the predictor 
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y£+i = yn + h ^7-i,k vkf n 
k=o 

i 

x n+i = x n + h y n +i + k 2 ^ ( a -i,k+i " 7-i,k+i)^ kf n 

k=o 

and the corrector 


yn+i = yn+i + h y_i,q vq+lf n+i 

x n+i = x n+l + h2 (°-i,q.+l " 7o,q+i) V* 1 f n+1 


(30) 


(31) 


Summed Form 


The predictor is q 

yn+i = An + k / 7-i ; k+i^ k ^n 

q. ' 

x J+i = B n + h2 ^-i,k+2V k f n 
k=o j 

The corrector is 

yn+i = yn+l + 1 f q+i^ +1 ^n+i 

» 

* , 1,2 _q+i,,* 

x n+l = x n+i + hff -i,q+2 v f n+i 

j 

and, of course, 

A n+i = A n + h f n+i 
B n+i = B n + hAn+i 


(32) 


(33) 


( 3*0 


In all these predictor- corrector algorithms, the calculation of the dif- 
ference table is facilitated by noting that, on defining 


e 


= fn+x - f 


* 

n+i 


the differences are 


V k f n +i = V^n+i + e , k = 1(1) q 

(see ref. 1, p. 1 96). The predicted differences, vAf**!, are obtained, of 
course, directly from the definition: 
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^n+i ~ 1 - f n 

^ + 1 fn+l = V^n+i “ V^n , k = l(l)q 
CONCLUDING REMARKS 


The present report displays a variety of algorithms for starting and 
continuing the numerical integration of systems of ordinary differential equa- 
tions. The exhaustive testing of these algorithms, for the purpose of compar- 
ing their effectiveness, would be an expensiVe process. Fortunately, a great 
deal of relevant experience has been obtained in recent years in computing 
installations throughout the world. In the present writer's opinion the best 
compromise between the conflicting desiderata of speed, accuracy, and program- 
ming compactness can be achieved by the following choice: 

(1) Use the fourth-order methods for first-order equations, sixth-order 
methods for second-order equations (q = 4, 6, respectively) . 

(2) Use the iterated starter, iterated eight times. 

(3) Use the summed form of the predictor-corrector algorithm, in back- 
ward difference form. 

(4) Carry four extra significant decimal digits, in floating-point form, 
to control round-off errors. 

The effectiveness of the summed form of the predict or -correct or algo- 
rithms has long been known to astronomers (ref. 3), and additional evidence is 
furnished by Henrici (ref. 1, pp. 336-339). The iterated starter is somewhat 
less efficient than the bootstrap version, but is far simpler to program and 
is much more modest in its storage requirements. 

The use of backward differences in the forward integration is preferable 
to the use of backward ordinates for two reasons: (l) the backward ordinate 

formula tends to add nearly equal quantities of alternating sign, whereas the 
backward difference formula adds monotonically decreasing quantities; and 
(2) the availability of the difference table makes error estimation and 
automatic adjustment of the interval size a straightforward procedure. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., April 19, 1965 
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APPENDIX 


CONVERSION OF ORDINAIES TO DIFFERENCES 


The calculation of a table of backward differences from a table of ordi- 
nates, and its extension in either direction when higher order differences are 
neglected, is a familiar procedure to computers working with pencil and paper. 
The programming of such procedures for automatic computers is less familiar, 
and the purpose of this appendix is to give some typical algorithms in FORTRAN 
format, the most widely used scientific programming language. 

The algorithms exhibited here are written as though x, y, f, A, and B 
were scalars. When they are vectors, the algorithms are easily generalized 
by the addition of a second subscript and suitable additional "DO LOOPS" over 
the range of the second subscript (the dimension of the vectors). 

Since the FORTRAN language does not permit subscripts nor indices in 
loops to assume nonpositive values, certain logical artificialities appear in 
these algorithms. Experienced programmers should have no difficulty in 
removing them for other, less restricted, programming languages. 


Backward Starter 

All the versions of the backward starter discussed in the main body of 
the report produce the backward ordinates, f _p, p = 0(l) q. Suppose these are 
placed in the array L: 

L(p + m) = f_ p 

where m is an arbitrary, positive integer. Then the nested DO LOOP 

_D0 k = 1, q 

l = q + 1 - k 
r DO 3=1 , l 

n = m + q - j 

— L L(n + l) = L(n) L(n + l) 

will yield vPf 0 in L(p + m), p = 0(l)q, q > 1. The structure of the algo- 
rithm is illustrated by the following diagram, in which q = 3 and m = 0 (in 
violation of the FORTRAN restriction!): 
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It is evident from the diagram that increasing q by unity adds one row to 
each column and adds an additional column on the right with one row. 


Forward Starter 

All the versions of the forward starter discussed in the main body of the 
report produce the forward ordinates, fp, p = 0(l)q* Suppose these are placed 
in the array L : 

L(p + m) = fp 

Then the nested DO LOOP 

k = 1, q. 

I = q_ + 1 - k 

J = 1, * 

n = m + q - j 

L(n + l) = L(n + l) - L(n) 
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■will yield V^fp in L(p + m), p = 0(l)q. The diagram illustrates the case 

q = 3) m = 0. 


Initial values 

k = 1, 2 = 3' 

k = 2, 2 = 2 

ii 

uo 

V* 

C-J 

II 

H 


j = 1, n = 2 

OJ 

II 

rf 

•s 

J— l 

II 

•ro 

1 

j = 1, n = 2 

L(3) = f 3 

L(3) = L(3) - L(2) 

L(3) = L(3) - L(2) 

L(3) = L(3) - L(2) 


= f 3 - f 2 

— v£*3 “ V^2 

= v*f 3 - V*f 2 


= Vfs 

= V2f3 

= v 3 ^ 


j = 2, n = 1 


j = 2, n = 1 


L(2) = L(2) - L(l) L(2) = L(2) - L(l) 


L(l) = fi 


= f 2 - fl 
= Vf2 

j = 3, n = 0 
L(l) = L(l) - L(0) 
= fl - fo 
= Wx 


= Vf 2 - yfi 

= ^f 2 


L(0) = f c 


To obtain backward differences of f 0 it is necessary to make the 
assumption that 7^- +1 f = 0. Then v^f is constant: 

V% = V^f 0 = V q f , p = l(l)q 

Then the additional nested DO LOOP 

r = q - 1 

DO k = 1, r 

2 = q - k 

r do j = l, 2 


n = m + q - j 

— L L(n) = L(n) - L(n + l) 

yields V-^fo i n L(p + m), p = 0(l)q, the desired result. The diagram 
illustrates the case q = b. m = 0. 
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TABLE I 


' 7 PA 



7 

27 

127 

247 

72O7 

14407 

604807 

12096O7 

p/k 

0 ~~ 

1 

2 

3 

4 " 

5 

6 

7 

-1 

1 

1 

5 

9 

251 

475 

19087 

36799 

0 

1 

-1 

-1 

-1 

-19 

-27 

-863 

-1375 

1 

1 

-3 

5 

1 

11 

11 

271 

351 

2 

1 

-5 

23 

-9 

-19 

-11 

-191 

-191 

3 

1 

-7 

53 

-55 

251 

27 

271 

191 

4 

1 

-9 

95 

-161 

1901 

-475 

-863 

-351 

5 

1 

-11 

149 

-351 

6731 

-4277 

19087 

1375 

6 

1 

-13 

215 

-649 

17261 

-17739 

198721 

-36799 

7 

1 

-15 

293 

-1079 

36731 

-52261 

9^3759 

-434241 


TABLE II.- a p ^ k 



a 

a 

12a 

12a 

240a 

240a 

60480a 

60480a 

3628800a 

p/k 

0 

1 

2 

3 

4 

5 

6 

7 

8 

-1 

1 

0 

1 

1 

19 

" 18 

4315 

4125 

237671 

0 

1 

-1 

1 

0 

-1 

-1 

-221 

-190 

-9829 

1 

1 

-2 

13 

-1 

-1 

0 

31 

31 

1571 

2 

1 

-3 

37 

-14 

19 

1 

31 

0 

-289 

3 

1 

-4 

73 

-51 

299 

-18 

-221 

-31 

-289 

4 

1 

-5 

121 

-124 

1319 

-317 

4315 

190 

1571 

5 

1 

-6 

181 

-245 

3799 

-I636 

84199 

-4125 

-9829 

6 

1 

-7 

253 

-426 

8699 

-5435 

496471 

-88324 

237671 

7 

1 

-8 

337 

-679 

17219 

-14134 

1866091 

-584795 

5537111 
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TABLE III.- 7c!,p,m 


CM 

O 

II 

a 1 

q. = 1: 127 

q. = 

2: 

24? 

p/m 

0 

0 

l 

0 

1 

2 

0 

-l 

-7 

l 

-15 

b 

-i 

1 

-3 

-13 

-5 

-25 

-12 

i 

2 


-7 

-23 

-23 

-28 

-9 

3 




-33 

4 

-55 



q 

= 3 : 7207 



a = 

4: 14407 


p/m 

0 

1 

2 

3 

0 

1 

2 

3 

4 

0 

-469 

177 

-87 

19 

-965 

462 

-336 

146 

-27 

1 

-739 

-393 

63 

-11 

-1467 

-830 

192 

-66 

11 

2 

-709 

-783 

-327 

19 

-1429 

-1522 

-720 

82 

-11 

3 

-739 

-633 

-897 

-251 

-1451 

-1374 

-1632 

-610 

27 

b 

-469 

-1743 

873 

-1901 

-1413 

-1586 

-1104 

-1902 

-475 

5 





-1915 

962 

-6336 

3646 

-4277 




q 

= 5: 604807 



p/m 

0 

1 

2 

3 

4 

5 

0 

-^1393 

23719 

-22742 

14762 

-5449 

863 

1 

-61343 

-36215 

10774 

-5482 

1817 

-271 

2 

-60209 

-62969 

-32150 

5354 

-1417 

191 

3 

-606 71 

-59063 

-65834 

-28330 

2489 

-271 

4 

-60209 

-62297 

-54998 

-71254 

-24256 

863 

5 

-61343 

-55031 

-75242 

-37738 

-84199 

-19087 

6 

-^1393 

-175865 

231274 

-456982 

248567 

-198721 





q = £ 

1: I2096O7 



p/m 

0 

1 

2 

3 

4 

5 

6 

0 

-84l6l 

55688 

-66109 

57024 

-31523 

9976 

-1375 

1 

-122335 

-74536 

26813 

-17984 

8899 

-2648 

351 

2 

-120609 

-124792 

-67165 

14528 

-5699 

1528 

-191 

3 

-121151 

-H9272 

-128803 

-6o48o 

7843 

-1688 

191 

4 

-120769 

-122488 

-115261 

-135488 

-53795 

3832 

-351 

5 

-121311 

-118312 

-129859 

-102976 

-147773 

-46424 

1375 

6 

-119585 

-130936 

-89437 

-177984 

-54851 

-176648 

-36799 

7 

-157759 

I38OO8 

-903715 

1198528 

-1465949 

717928 

-434241 


24 



TABLE IV.- a, 


q,p,m 


q = 0: 

12a 

q = 1: 12a 

q = 

2: 

240a 

p/m 

0 

0 

1 

0 

1 

2 

0 

1 

1 

0 

19 

2 

-i 

1 

13 

12 

1 

239 

22 

-i 

2 


23 

14 

479 

242 

19 

3 




739 

422 

299 



q = 

3: 240a 



q = 

4: 60480a 


p/m 

0 

1 

2 

3 

0 

1 

2 

3 

4 

0 

18 

5 

-4 

1 

4315 

2144 

-2334 

1136 

-221 

1 

239 

22 

-1 

0 

60259 

5420 

-66 

-124 

31 

2 

480 

239 

22 

-1 

120991 

60104 

5730 

-376 

31 

3 

721 

476 

245 

18 

181471 

120836 

6o4l4 

5420 

-221 

4 

942 

793 

368 

817 

241699 

182576 

118626 

62624 

4315 

5 





306715 

220124 

225726 

75476 

84199 




q 

= 5: 60480a 



p/m 

0 

1 

2 

3 

4 

5 

0 

4125 

3094 

-4234 

3036 

-1171 

190 

1 

60290 

5265 

244 

-434 

186 

-31 

2 

120991 

60104 

5730 

-376 

31 

0 

3 

l8l440 

120991 

60104 

5730 

-376 

31 

4 

241889 

181626 

120526 

60724 

5265 

-190 

5 

302590 

240749 

184476 

116726 

63574 

4125 

6 

358755 

327340 

178874 

266976 

54851 

88324 





q. = 6: 

3628800a 




p/m 

0 

1 

2 

3 

4 

5 

6 

0 

237671 

244614 

-401475 

378740 

-217695 

70374 

-9829 

1 

3618971 

306474 

38205 

-57460 

34725 

-11286 

1571 

2 

7259171 

3607974 

339465 

-16780 

-2475 

1734 

-289 

3 

10886111 

7261194 

3601905 

349580 

-26895 

3594 

-289 

4 

14514911 

10888134 

7255125 

3612020 

339465 

-20826 

1571 

5 

18145571 

14503914 

IO921125 

7200140 

3667005 

306474 

-9829 

6 

21762971 

18214374 

14297505 

11265140 

6856125 

3873414 

237671 

7 

25639271 

20099274 

23205465 

5979020 

19583625 

1865034 

5537111 


25 

























TABLE VI.- P 


q.,p,m 


q = 0: 

2P 

q = 1: 63 

q. = 

2: 

243 

p/m 

0 

0 

1 

0 

1 

2 

1 

1 

2 

1 

7 

6 

-1 

2 


4 

8 

16 

32 

0 

3 




27 

54 

27 



q. =. 

3 : 3603 



9 = . 

4: 14403 


p/m 

0 

1 

2 

3 

0 

1 

2 

3 

4 

1 

97 

114 

-39 

8 

367 

540 

-282 

116 

-21 

2 

224 

528 

-48 

16 

848 

2304 

-480 

256 

-48 

3 

351 

972 

243 

54 

1323 

4212 

486 

540 

-81 

4 

448 

1536 

384 

512 

1792 

6l44 

1538 

2048 

0 

5 





2375 

7500 

3750 

2500 

1875 




q. 

= 5 : 302403 



p/m 

0 

1 

2 

3 

4 

5 

1 

7386 

12945 

-9132 

5646 

-2046 

321 

2 

17040 

52224 

-17760 

13056 

-4848 

768 

3 

26568 

94527 

-1944 

23490 

-7776 

1215 

4 

36096 

136704 

16896 

58368 

-7680 

1536 

5 

45750 

178125 

37500 

93750 

18750 

4125 

6 

53136 

233280 

23328 

176256 

11664 

46656 





q = 6 : 

1209603 



p/m 

0 

l 

2 

3 

4 

5 

6 

1 

28549 

57750 

-51453 

42484 

-23109 

7254 

-995 

2 

65728 

223488 

-107520 

100864 

-55872 

17664 

-2432 

3 

102465 

400950 

-64881 

170100 

-88209 

27702 

-3807 

4 

139264 

577536 

-9216 

335872 

-107520 

36864 

-5120 

5 

176125 

753 750 

46875 

512500 

-28125 

57750 

-6875 

6 

212544 

933120 

93312 

705024 

46656 

186624 

0 

7 

257593 

1051638 

324135 

585844 

439383 

129654 

175273 
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"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof ” 

—National Aeronautics and Space Act of 1958 
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